Libraries

library(moments)
library(gvlma)
library(raster)
library(lmtest)
library(car)
library(leaps)
library(tidyverse)
library(leaflet)
library(ggplot2)
library(ggthemes)
library(tidyr)
library(broom)
library(sp)
library(sf)
library(vegan)
library(dplyr)
library(spData)
library(stats)
library(evaluate)
library(spatialEco)
library(knitr)
library(insol)
library(tmap)
library(geodiv)
library(GmAMisc)
library(viridis)
library(rgdal)
library(lidR)
library(lattice)
library(rasterVis)
library(RColorBrewer)
library(ggcorrplot)
library(cowplot)
library(googleway)
library(ggrepel)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthdata)

knitr::opts_chunk$set(cache=TRUE, warning=FALSE, message=FALSE)  

Read in Field Data

CanopyData <- read.csv("proj_data/CanopyData.csv")
PlotData <- read.csv("proj_data/PlotData.csv")
TreeData <- read.csv("proj_data/TreeData.csv")
TreeHeightsandAge <-read.csv("proj_data/Plot_tree_stats.csv")
DBHsAndDensity <-read.csv("proj_data/DBH_Density.csv") 

Bighorn National Forest Shape File

nf <- read_sf("proj_data/S_USA.AdministrativeForest/S_USA.AdministrativeForest.shp") %>%
  st_transform(.,crs=32613)
st_crs(nf)
## Coordinate Reference System:
##   User input: EPSG:32613 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 13N",
##     BASEGEOGCRS["WGS 84",
##         ENSEMBLE["World Geodetic System 1984 ensemble",
##             MEMBER["World Geodetic System 1984 (Transit)"],
##             MEMBER["World Geodetic System 1984 (G730)"],
##             MEMBER["World Geodetic System 1984 (G873)"],
##             MEMBER["World Geodetic System 1984 (G1150)"],
##             MEMBER["World Geodetic System 1984 (G1674)"],
##             MEMBER["World Geodetic System 1984 (G1762)"],
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ENSEMBLEACCURACY[2.0]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 13N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-105,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Between 108°W and 102°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],
##         BBOX[0,-108,84,-102]],
##     ID["EPSG",32613]]
bighorn <- dplyr::filter(nf,FORESTNAME=="Bighorn National Forest") 

Load DEM Data

DEM <- raster("proj_data/DEM.tif") 

st_crs(DEM)
## Coordinate Reference System:
##   User input: NAD83 (with axis order normalized for visualization) 
##   wkt:
## GEOGCRS["NAD83 (with axis order normalized for visualization)",
##     DATUM["North American Datum 1983",
##         ELLIPSOID["GRS 1980",6378137,298.257222101004,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

Load in Wyoming Landuse Landcover Data

#Read in the LULC Raster 
lulc=raster("proj_data/lulc.tif")

# LULC Names 
lulc_names =read_csv("proj_data/LULC.csv") %>%
  dplyr::select(class=ECOLSYS_LU...2, ID=Value)

# Convert to a Raster
lulc=as.factor(lulc)

# Update the RAT with a Left Join
levels(lulc)=left_join(levels(lulc)[[1]],lulc_names)
table(values(lulc))
## 
##      118      137      143      144      145      148      149      150 
##    22673   836200   663184  2557916   128707   119059 19613274  1427948 
##      151      155      158      184      191      192      194      311 
##  8817589  6583883     6633   443324      236       37     3863   120101 
##      312      314      315      317      323      325      398      425 
##     6797     9155  6992848    15082   478110   606140     5195   315854 
##      427      438      439      445      457      484      489      490 
##   290918    26349    39130    53143      133     2735     2243   337674 
##      491      495      529      535      537      545      546      547 
##  5339946   303973   114130    16366       36      131     2014      302 
##      549      552      554      557      558      561      567      568 
##  1796364    10377    12192      924    17332   253014     4524   124415 
##      573      579      581      582      583      584 
##   470583   192348   192184    56191    12897      402

Rename, Summarize, and Group Data

canopy <- CanopyData %>%
  group_by(Plot_Number) %>%
  dplyr::select(-Sunflecks) %>% 
  summarise_if(is.numeric, list(mean = function(x) mean(x, na.rm = TRUE), 
                                      sd = function(x) sd(x, na.rm = TRUE)))

plots <- PlotData %>%
 dplyr::group_by(Plot_Number) %>%
 dplyr::summarize(lat=mean(Lat), lon=mean(Lon)) 

trees<- TreeData %>%
  summarize(Plot_Number, Tree_ID, Species_Scientific_Name, DBH, Aprox_Tree_Age) %>%
   dplyr:: group_by(Plot_Number) 

Table of Species by Plot

tree1<-trees %>% 
  dplyr::select(-Tree_ID, -DBH) %>% 
  dplyr::group_by(Plot_Number,Species_Scientific_Name) %>% 
  dplyr::summarise(count=n()) %>%
  tidyr::spread(Species_Scientific_Name,count) %>% 
  dplyr:: mutate_all(function(x) ifelse(is.na(x),0,x))

treecounts <-table(trees$Species_Scientific_Name)
view(tree1)
view(treecounts)

Species Occurence Throughout Plots

This bar graph displays the total observations of each species across all plots. The Pinus contorta or Lodgepole Pine was observed the most frequently with more than 350 sightings. The Pinus flexilis, Limber Pine, and the Pinus ponderosa, the Ponderosa Pine each only had one sighting across the plots.

tree_counts_bargraph<-barplot(treecounts,
        main= "Species Occurence Throughout Plots", 
        horiz=TRUE, 
        xlab="Occurence", 
        ylab= "Species",
        col = c("#DD8317", "#3A9B44", "#47ACC2", "#EACF4F", "#EA5F4f", "#DD1794","#FFC300"), 
        names.arg=c("Juniperus occidentalis", "Picea engelmannii","Pinus contorta", "Pinus ponderosa", "Pinus flexilis", "Populus tremuloides", "Pseudotsuga menziesii"),
        cex.names=0.2, legend=TRUE)

Calculate the Diversity

diversity=data.frame(Plot_Number=tree1$Plot_Number,div=diversity(tree1))

Species Diversity

This bar graph shows the diversity within each plot. Plot 2 high the highest diversity at 1.2366849, and plot 60 had the lowest diversity of 0.0836497.

spatial_div<-ggplot(data= diversity, mapping=aes(x=Plot_Number, y=div, fill=div))+
  geom_bar(stat = "identity")+
  scale_fill_viridis_c(option="turbo")+
  labs(x = "Plot Number", y = "Diversity", title ="Species Diversity Across Plots")
spatial_div

Create a Merged Dataset

Merged <- merge(canopy, plots) %>%
  merge(DBHsAndDensity)%>%
 # merge(TreeHeightsandAge) %>%
  merge(diversity) %>% 
  st_as_sf(coords = c(x="lon", y="lat")) %>%
  st_set_crs(.,value=4326) %>%
  st_transform(.,crs=st_crs(DEM))

Merged_With_TreeHeightsandAge <- merge(canopy, plots) %>%
  merge(TreeHeightsandAge) %>%
  st_as_sf(coords = c(x="lon", y="lat")) %>%
  st_set_crs(.,value=4326) %>%
  st_transform(.,crs=st_crs(DEM))
#4326 would be fine 

tree_coords <- merge(trees, plots)

Calculate TPI, Slope angle, Slope aspect and Elevation

tpi10m <- raster("proj_data/tpi10m.tif")
tpi100m <- raster("proj_data/tpi100m.tif")
tpi1000m <- raster("proj_data/tpi1000m.tif")


TRI <- raster("proj_data/TRI.tif")
slope <- raster("proj_data/slope.tif")
aspect <- raster("proj_data/aspect.tif")

eastwest <- raster("proj_data/eastwest.tif")
northsouth <- raster("proj_data/northsouth.tif")

Create Landforms at 10 Meter Resolution

SD10m <- sd(tpi10m[],na.rm=T) 
SD10m
## [1] 0.4458888
landform_sd10m <- reclassify(tpi10m, matrix(c(-Inf, -SD10m, 1,
                                              -SD10m, -SD10m/2, 2,
                                              -SD10m/2, 0, 3,
                                              0, SD10m/2, 4,
                                              SD10m/2, SD10m, 5,
                                              SD10m, Inf, 6),
                                            ncol = 3, byrow = T),
                             right = T)

# Turn it into categorical raster
landform_sd10m <- as.factor(landform_sd10m) 

rat_sd10m <- levels(landform_sd10m)[[1]]

rat_sd10m[["landform_sd10m"]] <- c('Valley', 'Lower Slope', 
                                   'Flat Area','Middle Slope', 
                                   'Upper Slope', 'Ridge')

levels(landform_sd10m) <- rat_sd10m


# Plot the classification
tpi10mplot <-rasterVis::levelplot(landform_sd10m, col.regions = rev(brewer.pal(6,'RdYlBu')),
                       labels = rat_sd10m$landcover,
                       main = "10M TPI Landform Classification",
                       colorkey=list(labels=list(at=1:6, labels=rat_sd10m[["landform_sd10m"]])))

tpi10mplot

histtpi10m <-hist(tpi10m[])

Create Landforms at 100 Meter Resolution

SD100m <- sd(tpi100m[],na.rm=T)
SD100m 
## [1] 2.491715
# Make Landform Classifications 
#Morphologic class De Reu et al. 2013;  Weiss (2001)
landform_sd100m <- reclassify(tpi100m, matrix(c(-Inf, -SD100m, 1,
                                                -SD100m, -SD100m/2, 2,
                                                -SD100m/2, 0, 3,
                                                0, SD100m/2, 4,
                                                SD100m/2, SD100m, 5,
                                                SD100m, Inf, 6),
                                              ncol = 3, byrow = T),
                              right = T)

# Turn it into categorical raster
landform_sd100m <- as.factor(landform_sd100m) 

rat_sd100m <- levels(landform_sd100m)[[1]]

rat_sd100m[["landform_sd100m"]] <- c('Valley', 'Lower Slope', 
                                     'Flat Area','Middle Slope', 
                                     'Upper Slope', 'Ridge')

levels(landform_sd100m) <- rat_sd100m
#writeRaster(landform_sd100m, file="proj_data/landform_sd100m.grd", overwrite=TRUE)

# Plot the classification
tpi100mplot<- rasterVis::levelplot(landform_sd100m, col.regions = rev(brewer.pal(6,'RdYlBu')),
                        labels = rat_sd100m$landcover,
                        main = "100m TPI Landform Classification",
                        colorkey=list(labels=list(at=1:6, labels=rat_sd100m[["landform_sd100m"]])))


tpi100mplot

histtpi100m <- histtpi100<-hist(tpi100m[])

Create Landforms at 1000 Meter Resolution

SD1000m <- sd(tpi1000m[],na.rm=T)
SD1000m 
## [1] 29.7036
# Make landform classes
#Morphologic class De Reu et al. 2013;  Weiss (2001)
landform_sd1000m <- reclassify(tpi1000m, matrix(c(-Inf, -SD1000m, 1,
                                                  -SD1000m, -SD1000m/2, 2,
                                                  -SD1000m/2, 0, 3,
                                                  0, SD1000m/2, 4,
                                                  SD1000m/2, SD1000m, 5,
                                                  SD1000m, Inf, 6),
                                                ncol = 3, byrow = T),
                               right = T)

# Turn it into categorical raster
landform_sd1000m <- as.factor(landform_sd1000m) 

rat_sd1000m <- levels(landform_sd1000m)[[1]]

rat_sd1000m[["landform_sd1000m"]] <- c('Valley', 'Lower Slope', 
                                       'Flat Area','Middle Slope', 
                                       'Upper Slope', 'Ridge')

levels(landform_sd1000m) <- rat_sd1000m 
#writeRaster(landform_sd1000m, file="proj_data/landform_sd1000m.grd", overwrite= TRUE) #tif wont work for as.factor 

# Plot the classification
tpi1000mplot<- levelplot(landform_sd1000m, col.regions = rev(brewer.pal(6,'RdYlBu')),
                         labels = rat_sd1000m$landcover,
                         main = "1000m TPI Landform Classification",
                         colorkey=list(labels=list(at=1:6, labels=rat_sd1000m[["landform_sd1000m"]])))

tpi1000mplot

histtpi1000m<-hist(tpi1000m[])

Stack All Topographic Variables

all_topo <- stack(DEM,
              TRI,
              slope, aspect, 
              eastwest, northsouth,
              tpi10m, tpi100m, tpi1000m)
            
plot(all_topo)

Combine and Clean all Data at 10 Meter Resolution

envi <- raster::extract(all_topo, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, file="proj_data/envi.tif", overwrite=TRUE
                        ) %>%     
  scale() #this makes the regression coefficients comparable 


envi<-as.data.frame(envi)

bind<- bind_cols(Merged,envi) 


#LULC Introduction 
envi2 <- raster::extract(lulc, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean, file="proj_data/envi2.tif", overwrite=TRUE
                         ) %>%
  as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

bighornData <- bind_cols(bind,as.data.frame(envi2)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
    # lulc==~143~999,
    # lulc==144~ 999,
    # lulc==145~ 999,
    # lulc==315~ 999,
    # lulc==427~ 999,
    # lulc==490~ 999,
    # lulc==491~ 999,
    # lulc==568~ 999,
   TRUE~1,
  )))
#  dplyr::select(-"X")

Not Scaled- For Individual Maps!

enviNS <- raster::extract(all_topo, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean,  file="proj_data/enviNS.tif", overwrite=TRUE)

enviNS<-as.data.frame(enviNS)

bindNS<- bind_cols(Merged,enviNS) 

treesandbindNS<- merge(trees, bindNS)

#Add LULC
envi2NS <- raster::extract(lulc, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean, file="proj_data/envi2NS.tif", overwrite=TRUE
                         ) %>%
  as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

bighornDataNS <- bind_cols(bindNS,as.data.frame(envi2NS)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
    # lulc==~143~999,
    # lulc==144~ 999,
    # lulc==145~ 999,
    # lulc==315~ 999,
    # lulc==427~ 999,
    # lulc==490~ 999,
    # lulc==491~ 999,
    # lulc==568~ 999,
   TRUE~1,
  )))

Tree Height Combine and Clean at 10 Meter Resolution

enviTree <- raster::extract(all_topo, st_buffer(Merged_With_TreeHeightsandAge,
                                                dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, 
                            file="proj_data/enviTree.tif", overwrite=TRUE
                            )  %>%   
  scale() #this makes the regression coefficients comparable 
enviTree<-as.data.frame(enviTree) #%>% 


bindTree<- dplyr::bind_cols(Merged_With_TreeHeightsandAge,enviTree) %>%
  dplyr::rename(Average_Stand_Age=Average_Stamd_Age)

                

# LULC
envi2Tree <- raster::extract(lulc, st_buffer(Merged_With_TreeHeightsandAge, dist = 5, endCapStyle="SQUARE"), 
                             small=TRUE, fun=mean,file="proj_data/envi2Tree.tif", overwrite=TRUE) %>%
                              as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

bighornDataTreeHeightsAndAge <- bind_cols(bindTree,as.data.frame(envi2Tree)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  dplyr::select(-"X") %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
    TRUE~1,
  )))

Tree Height Not Scaled for Graphing

enviTreeNS <- raster::extract(all_topo, st_buffer(Merged_With_TreeHeightsandAge,
                                                dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, 
                            file="proj_data/enviTreeNS.tif", overwrite=TRUE
                            )  

enviTreeNS<-as.data.frame(enviTreeNS) 

bindTreeNS<- dplyr::bind_cols(Merged_With_TreeHeightsandAge,enviTreeNS) %>%
  dplyr::rename(Average_Stand_Age=Average_Stamd_Age)


#LULC
envi2TreeNS <- raster::extract(lulc, st_buffer(Merged_With_TreeHeightsandAge, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean,
                               file="proj_data/envi2TreeNS.tif", overwrite=TRUE) %>%
                              as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

NSbighornDataTreeHeightsAndAge <- bind_cols(bindTreeNS,as.data.frame(envi2TreeNS)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  dplyr::select(-"X") %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
    TRUE~1,
  )))

Forest Structure Correlations at 10 Meters

view(bighornData)

cor_bighorn10m <- bighornData %>% data.frame() %>%
  dplyr::select(
#"Plot_Number"      ,
"Gap_Fraction_mean",
"Stand_Density",
"Average_DBH",
"DBH_Sd",
"PAR_LAI_mean"     ,
"PAR_Average_mean"  ,
"Canopy_Density_mean",
"Leaf_Angle_mean" ,
"Gap_Fraction_sd",
"PAR_LAI_sd"  ,
    "PAR_Average_sd"    ,
    "Canopy_Density_sd",
    "Leaf_Angle_sd"   ,
    "div",
    "DEM"     ,
    "TRI"          ,
    "slope"        ,
    "aspect"  ,
    "eastwest"    ,
    "northsouth"       ,
    "tpi10m"          ,
    "tpi100m"          ,
    "tpi1000m"         , 
  )  %>%
  na.omit()

my_cor10<-cor(cor_bighorn10m)
view(my_cor10)
my_pmat10<-cor_pmat(cor_bighorn10m)
view(my_pmat10)
ggcorrplot(corr = my_cor10, p.mat = my_pmat10)

Forest Structure Regressions at 10 Meters

#Gap Fraction mean
#p-value: 2.378e-05
GapFraction_mean_10m<-lm(Gap_Fraction_mean~ TRI+ tpi1000m+lulc_reduced,
                         data=bighornData)
summary(GapFraction_mean_10m)
## 
## Call:
## lm(formula = Gap_Fraction_mean ~ TRI + tpi1000m + lulc_reduced, 
##     data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90669 -0.16654  0.00984  0.20381  0.72175 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.11244    0.09494  11.717  < 2e-16 ***
## TRI              0.07958    0.03974   2.003   0.0493 *  
## tpi1000m         0.08333    0.04003   2.082   0.0412 *  
## lulc_reduced149  0.49340    0.10851   4.547 2.39e-05 ***
## lulc_reduced151  0.60939    0.14439   4.220 7.61e-05 ***
## lulc_reduced155  0.16545    0.15041   1.100   0.2753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3257 on 66 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.349,  Adjusted R-squared:  0.2997 
## F-statistic: 7.076 on 5 and 66 DF,  p-value: 2.378e-05
#Gap Fraction sd
#p-value: 0.05405
GapFraction_sd_10m<-lm(Gap_Fraction_sd~  DEM,
                       data=bighornData)
summary(GapFraction_sd_10m)
## 
## Call:
## lm(formula = Gap_Fraction_sd ~ DEM, data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22355 -0.13059 -0.06727  0.07720  0.78374 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.29241    0.02247  13.015   <2e-16 ***
## DEM         -0.04431    0.02262  -1.959   0.0541 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.192 on 71 degrees of freedom
## Multiple R-squared:  0.05128,    Adjusted R-squared:  0.03791 
## F-statistic: 3.837 on 1 and 71 DF,  p-value: 0.05405
#####################################################################################
#PAR LAI mean
#p-value: 0.1319
PAR_LAI_mean_10m<-lm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,
                       data=bighornData)
summary(PAR_LAI_mean_10m)
## 
## Call:
## lm(formula = PAR_LAI_mean ~ DEM + tpi10m + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4342 -0.9292 -0.2874  0.9322  4.3386 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       2.1971     0.4383   5.013 4.17e-06 ***
## DEM               0.3347     0.1989   1.683   0.0970 .  
## tpi10m            0.3052     0.1788   1.707   0.0925 .  
## lulc_reduced149   0.4209     0.4970   0.847   0.4000    
## lulc_reduced151   0.5008     0.6919   0.724   0.4717    
## lulc_reduced155   1.3725     0.7065   1.942   0.0563 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.507 on 67 degrees of freedom
## Multiple R-squared:  0.1164, Adjusted R-squared:  0.0505 
## F-statistic: 1.766 on 5 and 67 DF,  p-value: 0.1319
#PAR LAI sd
#0.02881
PAR_LAI_sd_10m<-lm(PAR_LAI_sd~ eastwest,
                   data=bighornData)
summary(PAR_LAI_sd_10m)
## 
## Call:
## lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73833 -0.27252 -0.06329  0.20379  1.57294 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.64220    0.04844  13.259   <2e-16 ***
## eastwest    -0.10883    0.04877  -2.231   0.0288 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4138 on 71 degrees of freedom
## Multiple R-squared:  0.06553,    Adjusted R-squared:  0.05237 
## F-statistic: 4.979 on 1 and 71 DF,  p-value: 0.02881
#####################################################################################
#PAR Average mean
#p-value: 0.008254
PAR_Average_mean_10m<-lm(PAR_Average_mean~  lulc_reduced,
                         data=bighornData)
summary(PAR_Average_mean_10m)
## 
## Call:
## lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -410.28 -127.40  -67.51  107.00  777.92 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       420.48      61.84   6.799 3.06e-09 ***
## lulc_reduced149  -227.87      69.76  -3.266   0.0017 ** 
## lulc_reduced151  -286.57      94.46  -3.034   0.0034 ** 
## lulc_reduced155  -219.35      97.78  -2.243   0.0281 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 214.2 on 69 degrees of freedom
## Multiple R-squared:  0.1557, Adjusted R-squared:  0.119 
## F-statistic: 4.241 on 3 and 69 DF,  p-value: 0.008254
#PAR Average SD
#p-value: 0.01623
PAR_Average_sd10m<-lm(PAR_Average_sd~ TRI+slope+ northsouth,
                   data=bighornData)
summary(PAR_Average_sd10m)
## 
## Call:
## lm(formula = PAR_Average_sd ~ TRI + slope + northsouth, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -237.92 -126.03  -39.14   79.30  485.49 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   160.80      20.49   7.849 3.77e-11 ***
## TRI           129.52      82.37   1.572   0.1204    
## slope        -163.04      81.71  -1.995   0.0500 *  
## northsouth    -45.35      21.29  -2.130   0.0367 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 175 on 69 degrees of freedom
## Multiple R-squared:  0.1377, Adjusted R-squared:  0.1002 
## F-statistic: 3.672 on 3 and 69 DF,  p-value: 0.01623
#####################################################################################
#Canopy Density Mean 
#p-value: 8.287e-06
Canopy_Density_mean_10m<-lm(Canopy_Density_mean~  DEM+ lulc_reduced,
                      data=bighornData)
summary(Canopy_Density_mean_10m)
## 
## Call:
## lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.178  -2.206   1.187   4.191  22.762 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       55.305      2.925  18.908  < 2e-16 ***
## DEM                3.612      1.328   2.721  0.00826 ** 
## lulc_reduced149   14.206      3.317   4.283 5.93e-05 ***
## lulc_reduced151   12.912      4.605   2.804  0.00658 ** 
## lulc_reduced155    8.106      4.715   1.719  0.09010 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.06 on 68 degrees of freedom
## Multiple R-squared:  0.3421, Adjusted R-squared:  0.3034 
## F-statistic: 8.841 on 4 and 68 DF,  p-value: 8.287e-06
#Canopy Density SD
#5.175e-05
Canopy_Density_sd_10m<-lm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced,
                            data=bighornData)
summary(Canopy_Density_sd_10m)
## 
## Call:
## lm(formula = Canopy_Density_sd ~ DEM + tpi100m + lulc_reduced, 
##     data = bighornData)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.244 -2.638 -1.137  1.459 14.491 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       9.9574     1.4123   7.051 1.23e-09 ***
## DEM              -1.3150     0.6577  -1.999  0.04963 *  
## tpi100m          -1.7507     0.6103  -2.869  0.00551 ** 
## lulc_reduced149  -5.2374     1.5973  -3.279  0.00165 ** 
## lulc_reduced151  -5.5366     2.2326  -2.480  0.01566 *  
## lulc_reduced155  -4.4992     2.3043  -1.953  0.05506 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.842 on 67 degrees of freedom
## Multiple R-squared:  0.3281, Adjusted R-squared:  0.278 
## F-statistic: 6.544 on 5 and 67 DF,  p-value: 5.175e-05
#
####################################################################################
#Leaf Angle mean
#p-value: 0.5935
Leaf_Angle_mean_10<-lm(Leaf_Angle_mean~ lulc_reduced,
                          data=bighornData)
summary(Leaf_Angle_mean_10)
## 
## Call:
## lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.627  -8.346   1.694   9.114  35.162 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       48.727      4.390  11.100   <2e-16 ***
## lulc_reduced149    5.660      4.952   1.143    0.257    
## lulc_reduced151    1.049      6.706   0.156    0.876    
## lulc_reduced155    1.311      6.941   0.189    0.851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.21 on 69 degrees of freedom
## Multiple R-squared:  0.02696,    Adjusted R-squared:  -0.01534 
## F-statistic: 0.6374 on 3 and 69 DF,  p-value: 0.5935
#Leaf Angle sd
#p-value: 0.001394
Leaf_Angle_sd_10<-lm(Leaf_Angle_sd~ tpi1000m+lulc_reduced,
                       data=bighornData)
summary(Leaf_Angle_sd_10)
## 
## Call:
## lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.8374  -6.5164  -0.5118   7.8988  24.1936 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      24.7175     2.9241   8.453 3.67e-12 ***
## tpi1000m         -4.1758     1.2449  -3.354  0.00131 ** 
## lulc_reduced149  -6.5279     3.3152  -1.969  0.05308 .  
## lulc_reduced151  -0.9586     4.4813  -0.214  0.83126    
## lulc_reduced155  -3.1830     4.6759  -0.681  0.49839    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.13 on 67 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2296, Adjusted R-squared:  0.1836 
## F-statistic: 4.991 on 4 and 67 DF,  p-value: 0.001394
#####################################################################################
#Diversity
# p-value: 0.0009821
Diversity_10m<-lm(div~DEM+northsouth+lulc_reduced,
                     data=bighornData)
summary(Diversity_10m) 
## 
## Call:
## lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3811 -0.1322 -0.0067  0.1049  0.7027 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.38439    0.06148   6.252 3.22e-08 ***
## DEM              0.09229    0.02705   3.412   0.0011 ** 
## northsouth      -0.04444    0.02512  -1.769   0.0815 .  
## lulc_reduced149  0.16987    0.06971   2.437   0.0175 *  
## lulc_reduced151  0.03923    0.09604   0.408   0.6842    
## lulc_reduced155  0.15174    0.09824   1.545   0.1272    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.205 on 67 degrees of freedom
## Multiple R-squared:  0.2594, Adjusted R-squared:  0.2041 
## F-statistic: 4.694 on 5 and 67 DF,  p-value: 0.0009821
#####################################################################################

# Average DBH
#p-value: 0.0009843
DBH_mean_10m<-lm(Average_DBH~ tpi100m+ lulc_reduced,
                  data=bighornData)
summary(DBH_mean_10m) 
## 
## Call:
## lm(formula = Average_DBH ~ tpi100m + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.128  -4.770  -1.984   5.161  29.318 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       31.139      2.574  12.099  < 2e-16 ***
## tpi100m           -2.986      1.090  -2.739  0.00786 ** 
## lulc_reduced149   -6.592      2.900  -2.273  0.02616 *  
## lulc_reduced151    3.400      3.933   0.865  0.39026    
## lulc_reduced155   -6.044      4.172  -1.449  0.15202    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.904 on 68 degrees of freedom
## Multiple R-squared:  0.2353, Adjusted R-squared:  0.1903 
## F-statistic: 5.231 on 4 and 68 DF,  p-value: 0.0009843
#####################################################################################

#SD DBH
#p-value: 0.001216
DBH_sd_10m<-lm(DBH_Sd~ tpi1000m+ lulc_reduced,
                 data=bighornData)
summary(DBH_sd_10m) 
## 
## Call:
## lm(formula = DBH_Sd ~ tpi1000m + lulc_reduced, data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.6370  -2.5779  -0.6103   3.1628  19.7153 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       7.3101     1.5925   4.590    2e-05 ***
## tpi1000m         -1.2644     0.6780  -1.865  0.06658 .  
## lulc_reduced149   0.5346     1.8055   0.296  0.76805    
## lulc_reduced151   7.8764     2.4406   3.227  0.00194 ** 
## lulc_reduced155  -0.8562     2.5466  -0.336  0.73776    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.517 on 67 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.233,  Adjusted R-squared:  0.1872 
## F-statistic: 5.089 on 4 and 67 DF,  p-value: 0.001216
#####################################################################################
#Stand Density
#p-value: 1.806e-06
Stand_Density_10m<-lm(Stand_Density~ DEM+ tpi100m+ lulc_reduced,
               data=bighornData)
summary(Stand_Density_10m) 
## 
## Call:
## lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9979 -1.8267 -0.1134  1.6456  9.6533 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.0015     0.9799   6.125 5.38e-08 ***
## DEM               1.2977     0.4563   2.844  0.00591 ** 
## tpi100m           1.3722     0.4234   3.241  0.00186 ** 
## lulc_reduced149   3.0350     1.1082   2.739  0.00790 ** 
## lulc_reduced151  -0.2019     1.5491  -0.130  0.89671    
## lulc_reduced155   1.5211     1.5988   0.951  0.34483    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.36 on 67 degrees of freedom
## Multiple R-squared:  0.3969, Adjusted R-squared:  0.3519 
## F-statistic: 8.819 on 5 and 67 DF,  p-value: 1.806e-06
#####################################################################################

Tree Height Correlations at 10 Meeters

cor_bighornTreeHeightsAndAge <-bighornDataTreeHeightsAndAge %>% 
  data.frame() %>%
  dplyr::select(
"Mean_Tree_Height" ,
  "Sd_Tree_Height",
   "Max_Tree_Height",
    "Min_Tree_Height",
   "Average_Stand_Age",
    #"div",
    "DEM"     ,
    "TRI"          ,
    "slope"        ,
    "aspect"  ,
    "eastwest"    ,
    "northsouth"       ,
    "tpi10m"          ,
    "tpi100m"          ,
    "tpi1000m"         ,
  ) %>% 

  na.omit()

my_corTrees<-cor(cor_bighornTreeHeightsAndAge)
my_pmatTrees<-cor_pmat(cor_bighornTreeHeightsAndAge)
ggcorrplot(corr = my_corTrees,p.mat = my_pmatTrees)

Tree Height Regressions at 10 Meters

#Mean Tree Height
#p-value: 0.0199
 Mean_Tree_Height10m<-lm(Mean_Tree_Height~northsouth+lulc_reduced,
                     data=bighornDataTreeHeightsAndAge)
summary(Mean_Tree_Height10m)
## 
## Call:
## lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9490 -2.0172 -0.4666  2.3026  9.2168 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       7.3446     1.3570   5.412 3.64e-06 ***
## northsouth       -0.8025     0.5636  -1.424  0.16268    
## lulc_reduced149   5.3506     1.5301   3.497  0.00122 ** 
## lulc_reduced151   5.5459     2.3165   2.394  0.02170 *  
## lulc_reduced155   6.9510     3.6681   1.895  0.06572 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.43 on 38 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.259,  Adjusted R-squared:  0.181 
## F-statistic: 3.321 on 4 and 38 DF,  p-value: 0.0199
###################################################################################
#SD Tree Height
#p-value: 0.006275
SD_Tree_Height10m<-lm(Sd_Tree_Height~ tpi100m+lulc_reduced,
                      data=bighornDataTreeHeightsAndAge)
summary(SD_Tree_Height10m)
## 
## Call:
## lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0817 -1.2138 -0.1103  0.5723  4.4815 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.0266     0.6550   7.674 3.07e-09 ***
## tpi100m          -0.8031     0.2634  -3.049  0.00417 ** 
## lulc_reduced149  -1.0080     0.7271  -1.386  0.17371    
## lulc_reduced151   0.7790     1.0901   0.715  0.47922    
## lulc_reduced155  -0.8551     1.8526  -0.462  0.64704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.732 on 38 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.308,  Adjusted R-squared:  0.2352 
## F-statistic: 4.229 on 4 and 38 DF,  p-value: 0.006275
##################################################################################
#Max Tree Height (northsouth aspect)
#p-value: 0.00969
Max_Tree_Height10m<-lm(Max_Tree_Height~ aspect+ tpi1000m +lulc_reduced,
                       data=bighornDataTreeHeightsAndAge)
summary(Max_Tree_Height10m)
## 
## Call:
## lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6333 -2.8068  0.3041  1.9522  6.8380 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      15.8762     1.3137  12.085 2.07e-14 ***
## aspect            0.8112     0.5026   1.614   0.1150    
## tpi1000m         -1.1283     0.5897  -1.913   0.0635 .  
## lulc_reduced149   2.3294     1.4470   1.610   0.1159    
## lulc_reduced151   5.5133     2.1752   2.535   0.0156 *  
## lulc_reduced155   6.5971     3.6305   1.817   0.0773 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.388 on 37 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.326,  Adjusted R-squared:  0.235 
## F-statistic:  3.58 on 5 and 37 DF,  p-value: 0.00969
##################################################################################
#Average Stand Age 
#p-value: 0.0001305
Average_Stand_Age10m<-lm(Average_Stand_Age~aspect+ tpi1000m +lulc_reduced,
                         data=bighornDataTreeHeightsAndAge)

summary(Average_Stand_Age10m)
## 
## Call:
## lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.865 -10.901  -2.243   8.560  56.634 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       74.491      5.781  12.885  < 2e-16 ***
## aspect             5.211      2.422   2.151  0.03510 *  
## tpi1000m          -6.481      2.449  -2.646  0.01016 *  
## lulc_reduced149  -19.671      6.534  -3.010  0.00369 ** 
## lulc_reduced151    1.676      8.864   0.189  0.85066    
## lulc_reduced155  -18.959      9.432  -2.010  0.04852 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.89 on 66 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3113, Adjusted R-squared:  0.2591 
## F-statistic: 5.965 on 5 and 66 DF,  p-value: 0.0001305

Study Site Plots

# Interactive Map of All Plots 
tmap::tmap_mode(mode='view')
tm_shape(bindTree)+
  tm_dots(col='black')
#GGPlot of Study Site
studysiteplots<-ggplot(plots)+
  geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Longitude", y ="Latitude",
       title ="Study Site Plots", color ="Plot Number") + 
  theme(plot.title = element_text(hjust = 0.5))
studysiteplots

#Isolate Interactive Lidar Plots 
LidarplotsInt <-bindTree[-c(7,12,13, 23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43, 51,57, 58,59, 65,66, 69),]

#Interactive Map of Lidar Plots
tmap::tmap_mode(mode='view')
tm_shape(LidarplotsInt, )+
  tm_dots(col='blue')
#Isolate standard lidar study plots 
Lidarplots <-plots[-c(7,12,13, 23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43,
                              51,57, 58,59, 65,66, 69),]

#GGPlot of Lidar Plots 
lidarplotsgg<-ggplot(Lidarplots)+
  geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number))+
  scale_color_viridis_c(option="turbo")+
       labs(x = "Longitude", y ="Latitude",
       title ="Lidar Plots", color ="Plot Number") + 
  theme(plot.title = element_text(hjust = 0.5))
lidarplotsgg

Plot of Land Use Land Cover

lulc_colors=data.frame(class=levels(lulc)[[1]]$class)

lulc_colors$col=c("#EBA998","#FFFB00", "#4080BF", "#989F60", "#AF5063", "#3200FF", "#A289A8", "#BE20DF", "#971452", "#6250AF" ,  "#8EB880", "#624074", "#A0B232", "#B232A0", "#FF00E8","#F1197C", "#7C2B4B", "#7C2B72", "#A0D092", "#D095DE", "#70808F", "#10D4EF", "#AF50A7", "#EF103D","#4F30CF", "#FFC300", "#BCFF00", "#00BBFF", "#2080DF", "#DF8E20", "#9A8A86", "#CF3050", "#708F75", "#8F708D", "#95BCDE",  "#8F8F70", "#FFCD00", "#EB98E1", "#543E6C", "#DFBA20", "#FF9300", "#00FF43", "#B25932", "#B2327F", "#543B1B", "#1B543B", "#1B3D54", "#805C73", "#B62B85", "#F119B2", "#F18E19", "#591AAF", "#E0FF00", "#E0FF90")

lulcplot<-rasterVis::gplot(lulc)+
  geom_raster(aes(fill=as.factor(value)))+
  scale_fill_manual(labels=lulc_colors$class,
                    values=lulc_colors$col,
                    name="Landcover Type")+
  coord_equal()+
  theme(legend.position = "bottom", legend.key.size=unit(.5, "line"), legend.text=element_text(size=8))+
  guides(fill=guide_legend(ncol=3,byrow=TRUE))

lulcplot

Regression Results

regression_results <-read.csv("proj_data/Regression_Results.csv")
view(regression_results)


#Gap Fraction Mean 
par(mfrow=c(2,3))
plot(GapFraction_mean_10m, 1:6)

gvlma(GapFraction_mean_10m)
## 
## Call:
## lm(formula = Gap_Fraction_mean ~ TRI + tpi1000m + lulc_reduced, 
##     data = bighornData)
## 
## Coefficients:
##     (Intercept)              TRI         tpi1000m  lulc_reduced149  
##         1.11244          0.07958          0.08333          0.49340  
## lulc_reduced151  lulc_reduced155  
##         0.60939          0.16545  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = GapFraction_mean_10m) 
## 
##                     Value p-value                   Decision
## Global Stat        12.885 0.01185 Assumptions NOT satisfied!
## Skewness            1.457 0.22742    Assumptions acceptable.
## Kurtosis            1.715 0.19028    Assumptions acceptable.
## Link Function       4.881 0.02716 Assumptions NOT satisfied!
## Heteroscedasticity  4.832 0.02794 Assumptions NOT satisfied!
# Coefficients:
#     (Intercept)              TRI         tpi1000m  lulc_reduced149 lulc_reduced151  lulc_reduced155    
#         1.11244          0.07958          0.08333          0.49340          0.60939          0.16545  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = GapFraction_mean_10m) 
# 
#                     Value p-value                   Decision
# Global Stat        12.885 0.01185 Assumptions NOT satisfied!
# Skewness            1.457 0.22742    Assumptions acceptable.
# Kurtosis            1.715 0.19028    Assumptions acceptable.
# Link Function       4.881 0.02716 Assumptions NOT satisfied!
# Heteroscedasticity  4.832 0.02794 Assumptions NOT satisfied!

#Gap Fraction SD
par(mfrow=c(2,3))
plot(GapFraction_sd_10m, 1:6)

gvlma(GapFraction_sd_10m)
## 
## Call:
## lm(formula = Gap_Fraction_sd ~ DEM, data = bighornData)
## 
## Coefficients:
## (Intercept)          DEM  
##     0.29241     -0.04431  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = GapFraction_sd_10m) 
## 
##                      Value   p-value                   Decision
## Global Stat        75.8689 1.332e-15 Assumptions NOT satisfied!
## Skewness           37.5201 9.047e-10 Assumptions NOT satisfied!
## Kurtosis           36.6137 1.440e-09 Assumptions NOT satisfied!
## Link Function       0.5129 4.739e-01    Assumptions acceptable.
## Heteroscedasticity  1.2221 2.690e-01    Assumptions acceptable.
# Coefficients:
# (Intercept)          DEM  
#     0.29241     -0.04431  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = GapFraction_sd_10m) 
# 
#                      Value   p-value                   Decision
# Global Stat        75.8689 1.332e-15 Assumptions NOT satisfied!
# Skewness           37.5201 9.047e-10 Assumptions NOT satisfied!
# Kurtosis           36.6137 1.440e-09 Assumptions NOT satisfied!
# Link Function       0.5129 4.739e-01    Assumptions acceptable.
# Heteroscedasticity  1.2221 2.690e-01    Assumptions acceptable.


#PAR LAI Mean
par(mfrow=c(2,3))
plot(PAR_LAI_mean_10m, 1:6)

gvlma(PAR_LAI_mean_10m)
## 
## Call:
## lm(formula = PAR_LAI_mean ~ DEM + tpi10m + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM           tpi10m  lulc_reduced149  
##          2.1971           0.3347           0.3052           0.4209  
## lulc_reduced151  lulc_reduced155  
##          0.5008           1.3725  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = PAR_LAI_mean_10m) 
## 
##                      Value  p-value                   Decision
## Global Stat        9.45006 0.050784    Assumptions acceptable.
## Skewness           7.28542 0.006952 Assumptions NOT satisfied!
## Kurtosis           1.12472 0.288904    Assumptions acceptable.
## Link Function      1.01848 0.312880    Assumptions acceptable.
## Heteroscedasticity 0.02143 0.883608    Assumptions acceptable.
# Call:
# lm(formula = PAR_LAI_mean ~ DEM + tpi10m + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM           tpi10m  lulc_reduced149  lulc_reduced151  
#          2.1971           0.3347           0.3052           0.4209           0.5008  
# lulc_reduced155  
#          1.3725  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = PAR_LAI_mean_10m) 
# 
#                      Value  p-value                   Decision
# Global Stat        9.45006 0.050784    Assumptions acceptable.
# Skewness           7.28542 0.006952 Assumptions NOT satisfied!
# Kurtosis           1.12472 0.288904    Assumptions acceptable.
# Link Function      1.01848 0.312880    Assumptions acceptable.
# Heteroscedasticity 0.02143 0.883608    Assumptions acceptable.

#PAR LAI SD
par(mfrow=c(2,3))
plot(PAR_LAI_sd_10m, 1:6)

gvlma(PAR_LAI_sd_10m)
## 
## Call:
## lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
## 
## Coefficients:
## (Intercept)     eastwest  
##      0.6422      -0.1088  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = PAR_LAI_sd_10m) 
## 
##                      Value   p-value                   Decision
## Global Stat        31.2255 2.754e-06 Assumptions NOT satisfied!
## Skewness           15.6306 7.700e-05 Assumptions NOT satisfied!
## Kurtosis           14.2828 1.573e-04 Assumptions NOT satisfied!
## Link Function       0.7860 3.753e-01    Assumptions acceptable.
## Heteroscedasticity  0.5261 4.683e-01    Assumptions acceptable.
# Call:
# lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
# 
# Coefficients:
# (Intercept)     eastwest  
#      0.6422      -0.1088  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = PAR_LAI_sd_10m) 
# 
#                      Value   p-value                   Decision
# Global Stat        31.2255 2.754e-06 Assumptions NOT satisfied!
# Skewness           15.6306 7.700e-05 Assumptions NOT satisfied!
# Kurtosis           14.2828 1.573e-04 Assumptions NOT satisfied!
# Link Function       0.7860 3.753e-01    Assumptions acceptable.
# Heteroscedasticity  0.5261 4.683e-01    Assumptions acceptable.


#PAR Average Mean
par(mfrow=c(2,3))
plot(PAR_Average_mean_10m, 1:6)

gvlma(PAR_Average_mean_10m)
## 
## Call:
## lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
##           420.5           -227.9           -286.6           -219.4  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = PAR_Average_mean_10m) 
## 
##                        Value   p-value                   Decision
## Global Stat        4.265e+01 1.222e-08 Assumptions NOT satisfied!
## Skewness           1.811e+01 2.089e-05 Assumptions NOT satisfied!
## Kurtosis           1.321e+01 2.779e-04 Assumptions NOT satisfied!
## Link Function      1.853e-14 1.000e+00    Assumptions acceptable.
## Heteroscedasticity 1.133e+01 7.623e-04 Assumptions NOT satisfied!
# Call:
# lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#           420.5           -227.9           -286.6           -219.4  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = PAR_Average_mean_10m) 
# 
#                        Value   p-value                   Decision
# Global Stat        4.265e+01 1.222e-08 Assumptions NOT satisfied!
# Skewness           1.811e+01 2.089e-05 Assumptions NOT satisfied!
# Kurtosis           1.321e+01 2.779e-04 Assumptions NOT satisfied!
# Link Function      1.853e-14 1.000e+00    Assumptions acceptable.
# Heteroscedasticity 1.133e+01 7.623e-04 Assumptions NOT satisfied!

#PAR Average SD
par(mfrow=c(2,3))
plot(PAR_Average_sd10m, 1:6)

gvlma(PAR_Average_sd10m)
## 
## Call:
## lm(formula = PAR_Average_sd ~ TRI + slope + northsouth, data = bighornData)
## 
## Coefficients:
## (Intercept)          TRI        slope   northsouth  
##      160.80       129.52      -163.04       -45.35  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = PAR_Average_sd10m) 
## 
##                      Value   p-value                   Decision
## Global Stat        19.3889 0.0006590 Assumptions NOT satisfied!
## Skewness           12.4021 0.0004289 Assumptions NOT satisfied!
## Kurtosis            0.2763 0.5991333    Assumptions acceptable.
## Link Function       2.4977 0.1140115    Assumptions acceptable.
## Heteroscedasticity  4.2128 0.0401191 Assumptions NOT satisfied!
# Call:
# lm(formula = PAR_Average_sd ~ TRI + slope + northsouth, data = bighornData)
# 
# Coefficients:
# (Intercept)          TRI        slope   northsouth  
#      160.80       129.52      -163.04       -45.35  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = PAR_Average_sd10m) 
# 
#                      Value   p-value                   Decision
# Global Stat        19.3889 0.0006590 Assumptions NOT satisfied!
# Skewness           12.4021 0.0004289 Assumptions NOT satisfied!
# Kurtosis            0.2763 0.5991333    Assumptions acceptable.
# Link Function       2.4977 0.1140115    Assumptions acceptable.
# Heteroscedasticity  4.2128 0.0401191 Assumptions NOT satisfied!


#Canopy Density Mean
par(mfrow=c(2,3))
plot(Canopy_Density_mean_10m, 1:6)

gvlma(Canopy_Density_mean_10m)
## 
## Call:
## lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM  lulc_reduced149  lulc_reduced151  
##          55.305            3.612           14.206           12.912  
## lulc_reduced155  
##           8.106  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Canopy_Density_mean_10m) 
## 
##                     Value   p-value                   Decision
## Global Stat        87.465 0.000e+00 Assumptions NOT satisfied!
## Skewness           23.079 1.555e-06 Assumptions NOT satisfied!
## Kurtosis           51.359 7.694e-13 Assumptions NOT satisfied!
## Link Function       4.812 2.826e-02 Assumptions NOT satisfied!
## Heteroscedasticity  8.216 4.153e-03 Assumptions NOT satisfied!
# Call:
# lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#          55.305            3.612           14.206           12.912            8.106  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Canopy_Density_mean_10m) 
# 
#                     Value   p-value                   Decision
# Global Stat        87.465 0.000e+00 Assumptions NOT satisfied!
# Skewness           23.079 1.555e-06 Assumptions NOT satisfied!
# Kurtosis           51.359 7.694e-13 Assumptions NOT satisfied!
# Link Function       4.812 2.826e-02 Assumptions NOT satisfied!
# Heteroscedasticity  8.216 4.153e-03 Assumptions NOT satisfied!


#Canopy Density SD
par(mfrow=c(2,3))
plot(Canopy_Density_sd_10m, 1:6)

gvlma(Canopy_Density_sd_10m)
## 
## Call:
## lm(formula = Canopy_Density_sd ~ DEM + tpi100m + lulc_reduced, 
##     data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM          tpi100m  lulc_reduced149  
##           9.957           -1.315           -1.751           -5.237  
## lulc_reduced151  lulc_reduced155  
##          -5.537           -4.499  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Canopy_Density_sd_10m) 
## 
##                     Value   p-value                   Decision
## Global Stat        47.876 1.002e-09 Assumptions NOT satisfied!
## Skewness           20.297 6.629e-06 Assumptions NOT satisfied!
## Kurtosis            8.033 4.594e-03 Assumptions NOT satisfied!
## Link Function       3.087 7.890e-02    Assumptions acceptable.
## Heteroscedasticity 16.459 4.972e-05 Assumptions NOT satisfied!
# Call:
# lm(formula = Canopy_Density_sd ~ DEM + tpi100m + lulc_reduced, 
#     data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM          tpi100m  lulc_reduced149  lulc_reduced151  lulc_reduced155
#           9.957           -1.315           -1.751           -5.237           -5.537     -4.499  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Canopy_Density_sd_10m) 
# 
#                     Value   p-value                   Decision
# Global Stat        47.876 1.002e-09 Assumptions NOT satisfied!
# Skewness           20.297 6.629e-06 Assumptions NOT satisfied!
# Kurtosis            8.033 4.594e-03 Assumptions NOT satisfied!
# Link Function       3.087 7.890e-02    Assumptions acceptable.
# Heteroscedasticity 16.459 4.972e-05 Assumptions NOT satisfied!


#Leaf Angle Mean
par(mfrow=c(2,3))
plot(Leaf_Angle_mean_10, 1:6)

gvlma(Leaf_Angle_mean_10)
## 
## Call:
## lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
##          48.727            5.660            1.049            1.311  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Leaf_Angle_mean_10) 
## 
##                        Value p-value                Decision
## Global Stat        6.230e+00  0.1826 Assumptions acceptable.
## Skewness           2.686e+00  0.1012 Assumptions acceptable.
## Kurtosis           1.936e+00  0.1641 Assumptions acceptable.
## Link Function      5.268e-13  1.0000 Assumptions acceptable.
## Heteroscedasticity 1.608e+00  0.2048 Assumptions acceptable.
# Call:
# lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#          48.727            5.660            1.049            1.311  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Leaf_Angle_mean_10) 
# 
#                        Value p-value                Decision
# Global Stat        6.230e+00  0.1826 Assumptions acceptable.
# Skewness           2.686e+00  0.1012 Assumptions acceptable.
# Kurtosis           1.936e+00  0.1641 Assumptions acceptable.
# Link Function      5.268e-13  1.0000 Assumptions acceptable.
# Heteroscedasticity 1.608e+00  0.2048 Assumptions acceptable.

#Leaf Angle SD
par(mfrow=c(2,3))
plot(Leaf_Angle_sd_10, 1:6)

gvlma(Leaf_Angle_sd_10)
## 
## Call:
## lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)         tpi1000m  lulc_reduced149  lulc_reduced151  
##         24.7175          -4.1758          -6.5279          -0.9586  
## lulc_reduced155  
##         -3.1830  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Leaf_Angle_sd_10) 
## 
##                      Value p-value                Decision
## Global Stat        4.89145  0.2986 Assumptions acceptable.
## Skewness           0.48586  0.4858 Assumptions acceptable.
## Kurtosis           0.71337  0.3983 Assumptions acceptable.
## Link Function      3.67004  0.0554 Assumptions acceptable.
## Heteroscedasticity 0.02219  0.8816 Assumptions acceptable.
# Call:
# lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)         tpi1000m  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#         24.7175          -4.1758          -6.5279          -0.9586          -3.1830  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Leaf_Angle_sd_10) 
# 
#                      Value p-value                Decision
# Global Stat        4.89145  0.2986 Assumptions acceptable.
# Skewness           0.48586  0.4858 Assumptions acceptable.
# Kurtosis           0.71337  0.3983 Assumptions acceptable.
# Link Function      3.67004  0.0554 Assumptions acceptable.
# Heteroscedasticity 0.02219  0.8816 Assumptions acceptable

#Diversity
par(mfrow=c(2,3))
plot(Diversity_10m, 1:6)

gvlma(Diversity_10m)
## 
## Call:
## lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM       northsouth  lulc_reduced149  
##         0.38439          0.09229         -0.04444          0.16987  
## lulc_reduced151  lulc_reduced155  
##         0.03923          0.15174  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Diversity_10m) 
## 
##                      Value  p-value                   Decision
## Global Stat        17.0809 0.001864 Assumptions NOT satisfied!
## Skewness            6.2070 0.012725 Assumptions NOT satisfied!
## Kurtosis            5.4838 0.019194 Assumptions NOT satisfied!
## Link Function       0.8535 0.355558    Assumptions acceptable.
## Heteroscedasticity  4.5366 0.033178 Assumptions NOT satisfied!
# Call:
# lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM       northsouth  lulc_reduced149  lulc_reduced151  
#         0.38439          0.09229         -0.04444          0.16987          0.03923  
# lulc_reduced155  
#         0.15174  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Diversity_10m) 
# 
#                      Value  p-value                   Decision
# Global Stat        17.0809 0.001864 Assumptions NOT satisfied!
# Skewness            6.2070 0.012725 Assumptions NOT satisfied!
# Kurtosis            5.4838 0.019194 Assumptions NOT satisfied!
# Link Function       0.8535 0.355558    Assumptions acceptable.
# Heteroscedasticity  4.5366 0.033178 Assumptions NOT satisfied!

#DBH Mean
par(mfrow=c(2,3))
plot(DBH_mean_10m, 1:6)

gvlma(DBH_mean_10m)
## 
## Call:
## lm(formula = Average_DBH ~ tpi100m + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)          tpi100m  lulc_reduced149  lulc_reduced151  
##          31.139           -2.986           -6.592            3.400  
## lulc_reduced155  
##          -6.044  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = DBH_mean_10m) 
## 
##                     Value   p-value                   Decision
## Global Stat        22.860 0.0001351 Assumptions NOT satisfied!
## Skewness           10.835 0.0009959 Assumptions NOT satisfied!
## Kurtosis            8.499 0.0035529 Assumptions NOT satisfied!
## Link Function       1.769 0.1835180    Assumptions acceptable.
## Heteroscedasticity  1.757 0.1850515    Assumptions acceptable.
# Call:
# lm(formula = Average_DBH ~ tpi100m + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)          tpi100m  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#          31.139           -2.986           -6.592            3.400           -6.044  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = DBH_mean_10m) 
# 
#                     Value   p-value                   Decision
# Global Stat        22.860 0.0001351 Assumptions NOT satisfied!
# Skewness           10.835 0.0009959 Assumptions NOT satisfied!
# Kurtosis            8.499 0.0035529 Assumptions NOT satisfied!
# Link Function       1.769 0.1835180    Assumptions acceptable.
# Heteroscedasticity  1.757 0.1850515    Assumptions acceptable.

#DBH SD
par(mfrow=c(2,3))
plot(DBH_sd_10m, 1:6)

gvlma(DBH_sd_10m)
## 
## Call:
## lm(formula = DBH_Sd ~ tpi1000m + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)         tpi1000m  lulc_reduced149  lulc_reduced151  
##          7.3101          -1.2644           0.5346           7.8764  
## lulc_reduced155  
##         -0.8562  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = DBH_sd_10m) 
## 
##                      Value   p-value                   Decision
## Global Stat        23.8961 0.0000838 Assumptions NOT satisfied!
## Skewness            8.9941 0.0027085 Assumptions NOT satisfied!
## Kurtosis           10.2727 0.0013502 Assumptions NOT satisfied!
## Link Function       3.9400 0.0471502 Assumptions NOT satisfied!
## Heteroscedasticity  0.6893 0.4064005    Assumptions acceptable.
# Call:
# lm(formula = DBH_Sd ~ tpi1000m + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)         tpi1000m  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#          7.3101          -1.2644           0.5346           7.8764          -0.8562  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = DBH_sd_10m) 
# 
#                      Value   p-value                   Decision
# Global Stat        23.8961 0.0000838 Assumptions NOT satisfied!
# Skewness            8.9941 0.0027085 Assumptions NOT satisfied!
# Kurtosis           10.2727 0.0013502 Assumptions NOT satisfied!
# Link Function       3.9400 0.0471502 Assumptions NOT satisfied!
# Heteroscedasticity  0.6893 0.4064005    Assumptions acceptable.

#Stand Density
par(mfrow=c(2,3))
plot(Stand_Density_10m, 1:6)

gvlma(Stand_Density_10m)
## 
## Call:
## lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM          tpi100m  lulc_reduced149  
##          6.0015           1.2977           1.3722           3.0350  
## lulc_reduced151  lulc_reduced155  
##         -0.2019           1.5211  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Stand_Density_10m) 
## 
##                     Value p-value                Decision
## Global Stat        8.1492 0.08626 Assumptions acceptable.
## Skewness           2.6182 0.10564 Assumptions acceptable.
## Kurtosis           1.3827 0.23964 Assumptions acceptable.
## Link Function      3.3454 0.06739 Assumptions acceptable.
## Heteroscedasticity 0.8029 0.37023 Assumptions acceptable.
# Call:
# lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM          tpi100m  lulc_reduced149  lulc_reduced151  lulc_reduced155 
#          6.0015           1.2977           1.3722           3.0350          -0.2019   1.5211  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Stand_Density_10m) 
# 
#                     Value p-value                Decision
# Global Stat        8.1492 0.08626 Assumptions acceptable.
# Skewness           2.6182 0.10564 Assumptions acceptable.
# Kurtosis           1.3827 0.23964 Assumptions acceptable.
# Link Function      3.3454 0.06739 Assumptions acceptable.
# Heteroscedasticity 0.8029 0.37023 Assumptions acceptable.


# Tree Height Mean
par(mfrow=c(2,3))
plot(Mean_Tree_Height10m, 1:6)

gvlma(Mean_Tree_Height10m)
## 
## Call:
## lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)       northsouth  lulc_reduced149  lulc_reduced151  
##          7.3446          -0.8025           5.3506           5.5459  
## lulc_reduced155  
##          6.9510  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Mean_Tree_Height10m) 
## 
##                       Value p-value                Decision
## Global Stat        1.052881  0.9017 Assumptions acceptable.
## Skewness           0.997379  0.3179 Assumptions acceptable.
## Kurtosis           0.036870  0.8477 Assumptions acceptable.
## Link Function      0.005801  0.9393 Assumptions acceptable.
## Heteroscedasticity 0.012832  0.9098 Assumptions acceptable.
# Call:
# lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)       northsouth  lulc_reduced149  lulc_reduced151  lulc_reduced155 
#          7.3446          -0.8025           5.3506           5.5459          6.9510   
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Mean_Tree_Height10m) 
# 
#                       Value p-value                Decision
# Global Stat        1.052881  0.9017 Assumptions acceptable.
# Skewness           0.997379  0.3179 Assumptions acceptable.
# Kurtosis           0.036870  0.8477 Assumptions acceptable.
# Link Function      0.005801  0.9393 Assumptions acceptable.
# Heteroscedasticity 0.012832  0.9098 Assumptions acceptable.

# Tree Height SD
par(mfrow=c(2,3))
plot(SD_Tree_Height10m, 1:6)

gvlma(SD_Tree_Height10m)
## 
## Call:
## lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)          tpi100m  lulc_reduced149  lulc_reduced151  
##          5.0266          -0.8031          -1.0080           0.7790  
## lulc_reduced155  
##         -0.8551  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = SD_Tree_Height10m) 
## 
##                     Value p-value                   Decision
## Global Stat        4.9753 0.28984    Assumptions acceptable.
## Skewness           3.9299 0.04744 Assumptions NOT satisfied!
## Kurtosis           0.3955 0.52942    Assumptions acceptable.
## Link Function      0.3909 0.53182    Assumptions acceptable.
## Heteroscedasticity 0.2590 0.61080    Assumptions acceptable.
# Call:
# lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)          tpi100m  lulc_reduced149  lulc_reduced151  lulc_reduced155
#          5.0266          -0.8031          -1.0080           0.7790    -0.8551  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = SD_Tree_Height10m) 
# 
#                     Value p-value                   Decision
# Global Stat        4.9753 0.28984    Assumptions acceptable.
# Skewness           3.9299 0.04744 Assumptions NOT satisfied!
# Kurtosis           0.3955 0.52942    Assumptions acceptable.
# Link Function      0.3909 0.53182    Assumptions acceptable.
# Heteroscedasticity 0.2590 0.61080    Assumptions acceptable.

#Max Tree Height 
par(mfrow=c(2,3))
plot(Max_Tree_Height10m, 1:6)

gvlma(Max_Tree_Height10m)
## 
## Call:
## lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)           aspect         tpi1000m  lulc_reduced149  
##         15.8762           0.8112          -1.1283           2.3294  
## lulc_reduced151  lulc_reduced155  
##          5.5133           6.5971  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Max_Tree_Height10m) 
## 
##                      Value p-value                Decision
## Global Stat        4.26114 0.37182 Assumptions acceptable.
## Skewness           0.46840 0.49372 Assumptions acceptable.
## Kurtosis           0.85119 0.35622 Assumptions acceptable.
## Link Function      2.91962 0.08751 Assumptions acceptable.
## Heteroscedasticity 0.02193 0.88228 Assumptions acceptable.
# Call:
# lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced, 
#     data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)           aspect         tpi1000m  lulc_reduced149    lulc_reduced151  lulc_reduced155 
#         15.8762           0.8112          -1.1283           2.3294     5.5133           6.5971 
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Max_Tree_Height10m) 
# 
#                      Value p-value                Decision
# Global Stat        4.26114 0.37182 Assumptions acceptable.
# Skewness           0.46840 0.49372 Assumptions acceptable.
# Kurtosis           0.85119 0.35622 Assumptions acceptable.
# Link Function      2.91962 0.08751 Assumptions acceptable.
# Heteroscedasticity 0.02193 0.88228 Assumptions acceptable.

#Average Stand Age
par(mfrow=c(2,3))
plot(Average_Stand_Age10m, 1:6)

gvlma(Average_Stand_Age10m)
## 
## Call:
## lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)           aspect         tpi1000m  lulc_reduced149  
##          74.491            5.211           -6.481          -19.671  
## lulc_reduced151  lulc_reduced155  
##           1.676          -18.959  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Average_Stand_Age10m) 
## 
##                      Value  p-value                   Decision
## Global Stat        15.3846 0.003966 Assumptions NOT satisfied!
## Skewness            4.6874 0.030385 Assumptions NOT satisfied!
## Kurtosis            1.2152 0.270299    Assumptions acceptable.
## Link Function       0.8314 0.361869    Assumptions acceptable.
## Heteroscedasticity  8.6506 0.003270 Assumptions NOT satisfied!
# Call:
# lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced, 
#     data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)           aspect         tpi1000m  lulc_reduced149  lulc_reduced151  lulc_reduced155 
#          74.491            5.211           -6.481          -19.671         l1.676          -18.959     
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Average_Stand_Age10m) 
# 
#                      Value  p-value                   Decision
# Global Stat        15.3846 0.003966 Assumptions NOT satisfied!
# Skewness            4.6874 0.030385 Assumptions NOT satisfied!
# Kurtosis            1.2152 0.270299    Assumptions acceptable.
# Link Function       0.8314 0.361869    Assumptions acceptable.
# Heteroscedasticity  8.6506 0.003270 Assumptions NOT satisfied!

Gap Fraction Plots

#GapFraction_mean_10m<-lm(Gap_Fraction_mean~ TRI+ tpi1000m+lulc_reduced,data=bighornData)
GF.TRI<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=TRI, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Topographic Roughness Index", y ="Gap Fraction Mean",
       title ="Gap Fraction Mean vs Topographic Roughness IndexI", color ="Gap Fraction Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
GF.TRI

GF.100TPI<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "100m TPI", y ="Gap Fraction Mean",
       title ="Gap Fraction Mean vs 100m TPI", color ="Gap Fraction Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
GF.100TPI

GF.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Landcover Type", y ="Gap Fraction Mean",
       title ="Gap Fraction Mean vs Landcover Type", color ="Gap Fraction Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
GF.LULC

Gap_Fraction_mean_Car<-car::avPlots(model = GapFraction_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=TRUE,
             main=("Gap Fraction Added-Variable Plots")) 

#GapFraction_sd_10m<-lm(Gap_Fraction_sd~  DEM, data=bighornData)
GFsd.DEM<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Gap_Fraction_sd, color=Gap_Fraction_sd)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Elevation", y ="Gap Fraction Sd",
       title ="Gap Fraction Sd vs Elevation", color ="Gap Fraction Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
GFsd.DEM

Gap_Fraction_sd_Car<-car::avPlots(model = GapFraction_sd_10m, 
             col=carPalette()[1],
             col.lines=carPalette()[1], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Gap Fraction Added-Variable Plot")) 

PAR LAI Plots

#PAR LAI mean
#p-value: 0.1319
#PAR_LAI_mean_10m<-lm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,data=bighornData)

PARLAI.DEM<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "Elevation", y ="PAR Leaf Area Index Mean",
       title ="PAR Leaf Area Index Mean vs Elevation", color ="PAR LAI Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAI.DEM

PARLAI.TPI10<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi10m, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "10m TPI", y ="PAR Leaf Area Index Mean",
       title ="PAR Leaf Area Index Mean vs 10m TPI", color ="PAR LAI Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAI.TPI10

PARLAI.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "Landcover Type", y ="PAR Leaf Area Index Mean",
       title ="PAR Leaf Area Index Mean vs Landcover Type", color ="PAR LAI Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAI.LULC

PAR_LAI_CAR<-car::avPlots(model = PAR_LAI_mean_10m, 
             col=carPalette()[1],
             col.lines=carPalette()[1], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=TRUE,
             main=("PAR LAI Added-Variable Plots")) 

#PAR LAI sd
#p-value: 0.02881
#PAR_LAI_sd_10m<-lm(PAR_LAI_sd~ eastwest, data=bighornData)
PARLAIsd.EW<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=eastwest, y=PAR_LAI_sd, color=PAR_LAI_sd)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "East-West Aspect", y ="PAR Leaf Area Index Sd",
       title ="PAR Leaf Area Index Sd vs East-West Aspect", color ="PAR LAI Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAIsd.EW

PAR_LAI.sd_CAR<-car::avPlots(model = PAR_LAI_sd_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("PAR LAI Added-Variable Plot")) 

PAR Average Plots

#PAR Average mean
#p-value: 0.008254
#PAR_Average_mean_10m<-lm(PAR_Average_mean~  lulc_reduced,data=bighornData)

PARAvg.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=PAR_Average_mean, color=PAR_Average_mean)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Landcover Type", y ="PAR Average Mean",
       title ="PAR Average Mean vs Landcover Type", color ="PAR Average Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.LULC

PAR_Avg_CAR<-car::avPlots(model = PAR_Average_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("PAR Average mean Added-Variable Plots")) 

#PAR Average SD
#p-value: 0.01623
#PAR_Average_sd10m<-lm(PAR_Average_sd~ TRI+slope+ northsouth,data=bighornData)
PARAvg.sd.TRI<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=TRI, y=PAR_Average_sd, color=PAR_Average_sd)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Topographic Roughness Index", y ="PAR Average Sd",
       title ="PAR Average Sd vs Topographic Roughness Index", color ="PAR Average Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.TRI

PARAvg.sd.slope<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=slope, y=PAR_Average_sd, color=PAR_Average_sd)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Slope Anlge", y ="PAR Average Sd",
       title ="PAR Average Sd vs Slope Angle", color ="PAR Average Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.slope

PARAvg.sd.NS<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=northsouth, y=PAR_Average_sd, color=PAR_Average_sd)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "North-South Aspect", y ="PAR Average Sd",
       title ="PAR Average Sd vs North-South Aspect", color ="PAR Average Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.NS

PARAvg_Car<-car::avPlots(model = PAR_Average_sd10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("PAR Average sd Added-Variable Plots")) 

Canopy Density Plots

#Canopy Density Mean 
#p-value: 8.287e-06
#Canopy_Density_mean_10m<-lm(Canopy_Density_mean~  DEM+ lulc_reduced,data=bighornData)

CD.DEM <-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Canopy_Density_mean, color=Canopy_Density_mean)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Elevation", y ="Canopy Density Mean",
       title ="Canopy Density Mean vs Elevation", color ="Canopy Density Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
CD.DEM 

CD.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Canopy_Density_mean, color=Canopy_Density_mean)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Landcover Types", y ="Canopy Density Mean",
       title ="Canopy Density Mean vs Landcover Types", color ="Canopy Density Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
CD.LULC

CD_CAR<-car::avPlots(model = Canopy_Density_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Canopy Density mean Added-Variable Plots")) 

#Canopy Density SD
#5.175e-05
#Canopy_Density_sd_10m<-lm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced, data=bighornData)
CDs.DEM <-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Elevation", y ="Canopy Density Sd",
       title ="Canopy Density Sd vs Elevation", color ="Canopy Density Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
CDs.DEM 

CDs.TPI100<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "100m TPI", y ="Canopy Density Sd",
       title ="Canopy Density Sd vs 100m TPI", color ="Canopy Density Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
CDs.TPI100

CDs.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Landcover Types", y ="Canopy Density Sd",
       title ="Canopy Density Sd vs Landcover Types", color ="Canopy Density Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
CDs.LULC

CDs_CAR<-car::avPlots(model = Canopy_Density_sd_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Canopy Density sd Added-Variable Plots")) 

Leaf Angle Plots

#Leaf Angle mean
#p-value: 0.5935
#Leaf_Angle_mean_10<-lm(Leaf_Angle_mean~ lulc_reduced, data=bighornData)

LA.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Leaf_Angle_mean, color=Leaf_Angle_mean)) +
  scale_color_viridis_c(option="plasma")+
  labs(x = "Landcover Types", y ="Leaf Angle Mean",
       title ="Leaf Angle Mean vs Landcover Types", color ="Leaf Angle Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
LA.LULC

LA_CAR<-car::avPlots(model = Leaf_Angle_mean_10, 
             col=carPalette()[1],
             col.lines=carPalette()[1], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Leaf Angle mean Added-Variable Plots")) 

#Leaf Angle sd
#p-value: 0.001394
#Leaf_Angle_sd_10<-lm(Leaf_Angle_sd~ tpi1000m+lulc_reduced,data=bighornData)
LAs.TPI1000<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi1000m, y=Leaf_Angle_sd, color=Leaf_Angle_sd)) +
  scale_color_viridis_c(option="plasma")+
  labs(x = "1000m TPI", y ="Leaf Angle Sd",
       title ="Leaf Angle Sd vs 1000m TPI", color ="Leaf Angle Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
LAs.TPI1000

LAs.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Leaf_Angle_sd, color=Leaf_Angle_sd)) +
  scale_color_viridis_c(option="plasma")+
  labs(x = "Land Cover Types", y ="Leaf Angle Sd",
       title ="Leaf Angle Sd vs Land Cover Types", color ="Leaf Angle Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
LAs.LULC

LA.sd_CAR<-car::avPlots(model = Leaf_Angle_sd_10, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Leaf Angle sd Added-Variable Plots")) 

Diversity Plots

#Diversity
# p-value: 0.0009821
#Diversity_10m<-lm(div~DEM+northsouth+lulc_reduced, data=bighornData)

div.DEM<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=div, color=div)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Elevation", y ="Diversity",
       title ="Diversity vs Elevation", color ="Diversity") + 
  theme(plot.title = element_text(hjust = 0.5))
div.DEM

div.NS<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=northsouth, y=div, color=div)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "North-South Aspect", y ="Diversity",
       title ="Diversity vs North-South Aspect", color ="Diversity") + 
  theme(plot.title = element_text(hjust = 0.5))
div.NS

div.lulc<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=div, color=div)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Landcover Types", y ="Diversity",
       title ="Diversity vs Landcover Types", color ="Diversity") + 
  theme(plot.title = element_text(hjust = 0.5))
div.lulc

div_Car<-car::avPlots(model = Diversity_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Diversity Added-Variable Plots")) 

DBH Plots

# Average DBH
#p-value: 0.0009843
#DBH_mean_10m<-lm(Average_DBH~ tpi100m+ lulc_reduced,data=bighornData)
DBH_Avg.TPI100<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Average_DBH, color=Average_DBH)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "100m TPI", y ="Average DBH",
       title ="Average Diameter Breast Height vs 100m TPI", color ="Average DBH") + 
  theme(plot.title = element_text(hjust = 0.5))
DBH_Avg.TPI100

DBH_Avg.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Average_DBH, color=Average_DBH)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "Landcover Types", y ="Average DBH",
       title ="Average Diameter Breast Height vs Landcover Types", color ="Average DBH") + 
  theme(plot.title = element_text(hjust = 0.5))
DBH_Avg.LULC

DBH_Avg_CAR<-car::avPlots(model = DBH_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("DBH mean Added-Variable Plots")) 

#SD DBH
#p-value: 0.001216
#DBH_sd_10m<-lm(DBH_Sd~ tpi1000m+ lulc_reduced, data=bighornData)
dbh.sd.tpi1000<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi1000m, y=DBH_Sd, color=DBH_Sd)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "1000m TPI", y ="DBH Sd",
       title ="Sd of Diameter Breast Height vs 1000m TPI", color ="DBH Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
dbh.sd.tpi1000

dbh.sd.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=DBH_Sd, color=DBH_Sd)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "Land Cover Types", y ="DBH Sd",
       title ="Sd of Diameter Breast Height vs Land Cover Types", color ="DBH Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
dbh.sd.LULC

dbh.sd.CAR<-car::avPlots(model = DBH_sd_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("DBH sd Added-Variable Plots")) 

Stand Density Plots

#Stand Density
#p-value: 1.806e-06
#Stand_Density_10m<-lm(Stand_Density~ DEM+ tpi100m+ lulc_reduced, data=bighornData)

den.dem<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Stand_Density, color=Stand_Density)) +
  scale_color_viridis_c(option="mako")+
  labs(x = "Elevation", y ="Stand Density",
       title ="Stand Density vs Elevation", color ="Stand Density") + 
  theme(plot.title = element_text(hjust = 0.5))
den.dem

den.tpi100<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Stand_Density, color=Stand_Density)) +
  scale_color_viridis_c(option="mako")+
  labs(x = "100m TPI", y ="Stand Density",
       title ="Average Tree Height vs 100m TPI", color ="Stand Density") + 
  theme(plot.title = element_text(hjust = 0.5))
den.tpi100

den.lulc<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Stand_Density, color=Stand_Density)) +
  scale_color_viridis_c(option="mako")+
  labs(x = "Landcover Types", y ="Stand Density",
       title ="Stand Density vs Landcover Types", color ="Stand Density") + 
  theme(plot.title = element_text(hjust = 0.5))
den.lulc

den.car<-car::avPlots(model = Stand_Density_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Stand Density Added-Variable Plots")) 

Tree Height Plots

#Mean Tree Height
#p-value: 0.0199
#Mean_Tree_Height10m<-lm(Mean_Tree_Height~northsouth+lulc_reduced, data=bighornDataTreeHeightsAndAge)

avgTH.NS<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=northsouth, y=Mean_Tree_Height, color=Mean_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "North-South Aspect", y ="Average Tree Height",
       title ="Average Tree Height vs North-South Aspect", color ="Average Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
avgTH.NS

avgTH.LULC<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Mean_Tree_Height, color=Mean_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Land Cover Types", y ="Average Tree Height",
       title ="Average Tree Height vs Land Cover Types", color ="Average Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
avgTH.LULC

avgTH.CAR<-car::avPlots(model = Mean_Tree_Height10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Average Tree Height Added-Variable Plots")) 

#SD Tree Height
#p-value: 0.006275
#SD_Tree_Height10m<-lm(Sd_Tree_Height~ tpi100m+lulc_reduced,  data=bighornDataTreeHeightsAndAge)

sd.th.tpi100<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=tpi100m, y=Sd_Tree_Height, color=Sd_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "100m TPI", y ="Tree Height Sd",
       title ="Tree Height Sd vs 100m TPI", color ="Tree Height Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
sd.th.tpi100

sd.th.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Sd_Tree_Height, color=Sd_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Land Cover Types", y ="Tree Height Sd",
       title ="Tree Height Sd vs Aspect", color ="Tree Height Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
sd.th.lulc

sd.thCAR<-car::avPlots(model = SD_Tree_Height10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Tree Height sd Added-Variable Plots")) 

#Max Tree Height (northsouth aspect)
#p-value: 0.00969
#Max_Tree_Height10m<-lm(Max_Tree_Height~ aspect+ tpi1000m+lulc_reduced,data=bighornDataTreeHeightsAndAge)
max.th.aspect<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=aspect, y=Max_Tree_Height, color=Max_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Aspect", y ="Max Tree Height",
       title ="Maximum Tree Heights vs Aspect", color ="Maximum Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
max.th.aspect

max.th.tpi1000<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=tpi1000m, y=Max_Tree_Height, color=Max_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "1000m TPI", y ="Max Tree Height",
       title ="Maximum Tree Heights vs 1000m TPI", color ="Maximum Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
max.th.tpi1000

max.th.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Max_Tree_Height, color=Max_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Landcover Types", y ="Max Tree Height",
       title ="Maximum Tree Heights vs Landcover Type", color ="Maximum Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
max.th.lulc

max.th.CAR<-car::avPlots(model = Max_Tree_Height10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Maximum Tree Height Added-Variable Plots")) 

Stand Age Plots

#Average Stand Age 
#p-value: 0.0001305
#Average_Stand_Age10m<-lm(Average_Stand_Age~aspect+ tpi1000m +lulc_reduced,  data=bighornDataTreeHeightsAndAge)
avg.sa.aspect<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=aspect, y=Average_Stand_Age, color=Average_Stand_Age)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "aspect", y ="Average_Stand_Age", title ="Average Stand Age vs Aspect", color ="Average_Stand_Age") +
  theme(plot.title = element_text(hjust = 0.5))
avg.sa.aspect

avg.sa.tpi1000<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=tpi1000m, y=Average_Stand_Age, color=Average_Stand_Age)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "1000m TPI", y ="Average Stand Age", title ="Average Stand Age vs 1000m TPI", color ="Average Stand Age") +
  theme(plot.title = element_text(hjust = 0.5))
avg.sa.tpi1000

avg.sa.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Average_Stand_Age, color=Average_Stand_Age)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Landcover Types", y ="Average Stand Age",
       title ="Average Stand Age vs Landcover Type", color ="Average Stand Age") + 
  theme(plot.title = element_text(hjust = 0.5))
avg.sa.lulc

avg.sa.CAR<-car::avPlots(model = Average_Stand_Age10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Estimated Stand Age Added-Variable Plots")) 

Lidar Images of Plots

#### Read in Plots
plotsf<- plots %>%
  st_as_sf(coords = c(x="lon", y="lat")) %>%
  st_set_crs(.,value=4326) 

# Obtain Plot Coordinates in Lidar Projection
lasfile="proj_data/LAS/plot_52_USGS_LPC_WY_Sheridan_2020_D20_13TCK1344.laz"

lproj=readLASheader(lasfile) %>% 
  crs()

#10m polys 
plot_polys <- plotsf %>% 
  st_transform(.,crs="+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m no_defs") %>% 
  st_buffer(dist = 50,endCapStyle="SQUARE") %>% 
  st_transform(.,crs=st_crs(4326)) 

# Plot Specific Change 
plotgeo <-plot_polys[52,] %>% 
  st_transform(.,crs=lproj) 


# Read in New LAS File
lasPlots <-readLAS(lasfile,
                   filter = paste("-keep_xy ",paste(st_bbox(plotgeo),
                                                    collapse=" "))) %>% 
  clip_roi(geometry= plotgeo) 

plot(lasPlots) 

Plot 2.1, Lidar Image of Plot 2 in Bighorn National Forest. Plot 2.2, Lidar Image of Plot 2 in Bighorn National Forest. Plot 6, Lidar Image of Plot 6 in Bighorn National Forest. Plot 8.1, Lidar Image of Plot 8 in Bighorn National Forest.

Plot 8.2, Lidar Image of Plot 8 in Bighorn National Forest. Plot 11.1, Lidar Image of Plot 11 in Bighorn National Forest. Plot 11.2, Lidar Image of Plot 11 in Bighorn National Forest. Plot 14.1, Lidar Image of Plot 14 in Bighorn National Forest. Plot 14.2, Lidar Image of Plot 14 in Bighorn National Forest. Plot 15, Lidar Image of Plot 15 in Bighorn National Forest. Plot 17, Lidar Image of Plot 17 in Bighorn National Forest. Plot 18, Lidar Image of Plot 18 in Bighorn National Forest.

Plot 20, Lidar Image of Plot 20 in Bighorn National Forest.

Plot 42.1, Lidar Image of Plot 42 in Bighorn National Forest. Plot 42.2, Lidar Image of Plot 42 in Bighorn National Forest. Plot 46, Lidar Image of Plot 46 in Bighorn National Forest.

Plot 47.1, Lidar Image of Plot 47 in Bighorn National Forest. Plot 47.2, Lidar Image of Plot 47 at 30m Resolution in Bighorn National Forest. Plot 48.1, Lidar Image of Plot 48 at 30m Resolution in Bighorn National Forest. Plot 48.2, Lidar Image of Plot 48 at 30m Resolution in Bighorn National Forest. Plot 50, Lidar Image of Plot 50 in Bighorn National Forest. Plot 52.1, Lidar Image of Plot 52 at 100m Resolution in Bighorn National Forest. Plot 52.2, Lidar Image of Plot 52 at 30m Resolution in Bighorn National Forest.

Calculating Tree Heights

#dtm, dsm, chm 
dsm2 <- lidR::grid_canopy(lasPlots, res=2, p2r()) 
dtm2 <- lidR::grid_terrain(lasPlots, res= 2, algorithm=tin()) 

#Crop! 
chm2<-dsm2-dtm2 
chmcrop <- crop(chm2, plotgeo)  

#TRUE CROPPED MEAN 
chmmean<- cellStats(chmcrop, stat=mean) 
view(chmmean)  

#TRUE CROPPED SD
chmstd<-  cellStats(chmcrop, stat=sd) 
view(chmstd) 

#Max Value
chmmax <- cellStats(chmcrop, stat=max)
view(chmmax)

#Min Value 
chmmin <- cellStats(chmcrop, stat=min)
view(chmmin)

#Plot and Hist
plot(chmcrop, col = height.colors(25))

hist(chmcrop)

Plotted Species Diversity

Plotted_div<-ggplot(data=Merged, mapping=aes(x=plots$lon, y=plots$lat, color =div))+
  geom_point()+
  scale_color_viridis_c(option="turbo")+
   labs(x = "Longitude", y = "Latitude", title ="Species Occurene Throughout Plots", color ="Diversity")

Plotted_div

Statistics Per Plot

This image displays the data collected at each plot in a spatial format. The graphic in the top left shows the plots, numbered. The next image shows the mean Gap Fraction per plot. Followed by the mean PAR LAI per plot, the mean PAR Average per plot and the mean Canopy Density per plot. The second row of images shows the mean Leaf Angle per plot, the standard deviation of the Gap Fraction per plot, the standard deviation of the PAR LAI per plot, the standard deviation of the PAR Average per plot and the standard deviation of the Canopy Density per plot.

plot(bind) 
 title(main = "Extracting Points", line= 8)

Trees

This graph displays the location of species across plots, with a size of the point correlated to the diameter breast height of the individual trees.

spatial_tc<-ggplot(data=tree_coords, mapping=aes(x=lon, y=lat, color=Species_Scientific_Name, size=bighornDataNS$st))+
  geom_point()+  
  scale_color_viridis_d(option="turbo")+
  labs(x = "Latitude", y = "Longitude", title ="Spatial Occurrence of Tree Species", color ="Species Scientific Name") +
  theme(plot.title = element_text(hjust = 0.5))

spatial_tc